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Abstract 

We present a technique for the rapid and reliable prediction of linear-functional outputs of ellip- 
tic (and parabolic) partial differential equations with affine parameter dependence. The essential 
components are (i) (provably) rapidly convergent global reduced- basis approximations Galcrkin 
projection onto a space VVW spanned by solutions of the governing partial differential equation 
at N selected points in parameter space; («) a posteriori error estimation relaxations of the 
error-residual equation that provide inexpensive yet sharp and rigorous bounds for the error in the 
outputs of interest; and (Hi) off-line/on line computational procedures — methods which decouple 
the generation and projection stages of the approximation process. The operation count for the 
on-line stage — in which, given a new parameter value, we calculate the output of interest and 
associated error bound — ■ depends only on N (typically very small) and the parametric complexity 
of the problem; the method is thus ideally suited for the repeated and rapid evaluations required 
in the context of parameter estimation, design, optimization, and real-time control. 


1 Introduction 

The optimization, control, and characterization of an engineering component or system requites the 
prediction of certain “quantities of interest,” or performance metrics, which we shall denote outputs - 
for example deflections, maximum stresses, maximum temperatures, heat transfer rates, flowrates, or 
lift and drags. These outputs are typically expressed as functionals of field variables associated with 
a parametrized partial differential equation which describes the physical behavior of the component or 
system. The parameters, which we shall denote inputs , serve to identify a particular “configuration” 
of the component: these inputs may represent design or decision variables, such as geometry for 
example, in optimization studies; control variables, such as actuator power for example in real- 
time applications; or characterization variables, such as physical properties — for example in inverse 
problems. We thus arrive at an implicit input-output relationship, evaluation of which demands solution 
of the underlying partial differential equation. 

Our goal is the development of computational methods that permit rapid and reliable evaluation 
of this partial- differential-equation-induced input-output relationship in the limit of many queries 
that is, in the design, optimization, control, and characterization contexts. The “many query” limit has 
certainly received considerable attention: from “fast loads” or multiple right-hand side notions (e.g., [7, 
91) to matrix perturbation theories (e.g., [1, 28]) to continuation methods (e.g., [2, 23]). Our particular 
approach is based on the reduced-basis method, first introduced in the late 1970s for nonlinear structural 
analysis (3, 19], and subsequently developed more broadly in the 1980s and 1990s [5, 6, 10, 21, 22, 24]. 
The reduced-basis method recognizes that the field variable is not, in fact, some arbitrary member 
of the infinite-dimensional solution space associated with the partial differential equation; rather, it 
resides, or “evolves,” on a much lower-dimensional manifold induced by the parametric dependence. 

The reduced-basis approach as earlier articulated is local in parameter space in both practice and 
theory. To wit, Lagrangian or Taylor approximation spaces for the low-dimensional manifold are 
typically defined relative to a particular parameter point; and the associated a priori convergence theory 
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relies on asymptotic arguments m sufficiently small neighborhoods (10). As a result, the computational 
improvements relative to conventional (say) finite element approximation ~~ are often quite modest 
[22]. Our work differs from these earlier efforts in several important ways: first, we develop (in some 
cases, provably) global approximation spaces; second, we introduce rigorous a posteriori error estimators’ 
and third, we exploit off-line/ on-line computational decompositions (see [5j for an earlier application 
ol this stiategy within the reduced basis context). These three ingredients allow us — for the restricted 
but important class of “parameter-affine” problems — to reliably decouple the generation and projection 
stages of reduced -basis approximation, thereby effecting computational economies of several orders of 
magnitude. 

In this expository review paper wc focus on these new ingredients. In Section 2 we introduce an 
abstract problem formulation and several illustrative instantiations. In Section 3 we describe, for coer- 
cive symmetric problems and “compliant” outputs, the reduced -basis approximation; and in Section 4 
we present the associated a posteriori error estimation procedures. In Section 5 we consider the exten- 
sion of our approach to noncompliant outputs and nonsymmetric operators; eigenvalue problems* and 
more briefly, noncoercive operators, parabolic equations, and non-affine problems. A description of the 
system architecture iu which these numerical objects reside may be found in [26 j . 


2 Problem Statement 


2.1 Abstract Formulation 


We consider a suitably regular domain 0 C It d , d = 1, 2, or 3, and associated function space X C 
where H 1 ^) = {u € Vv e (L 2 (0))‘ i }, and L 2 (H) is the space of square integrable 

functions over fh The inner product and norm associated with X are given by (*, -) y and |j*jj \ = (• .)*/2 
respectively. We also define a parameter set D 6 IR^, a particular point in which will be denoted p. 
Note that Q does not depend on the parameter. 

We then introduce a bilinear form a: X x X x V 11, and linear forms /: X -> JR, £ : X -> IR. 
We shall assume that a is continuous, a(w } v\p) < 7 (/r) H|x IMU < 7o \\w\\ x \\v\\ Xl Vp 6 v\ 
fur thei more, in Sections 3 and 4, we assume that u is coercive, 


0 < ao < 


a(/r) = inf 
wjG X 


a(tn, w ; p) 

IMIx 


V/ieD, 


(i) 


and symmetric, a('WyV\p) — a(v y w;p) } Vu;,v eI,V^D. We also require that our linear forms / and 
£ be bounded; in Sections 3 and 4 we additionally assume a “compliant” output, f(v) = £(v)y V?; G X 
We shall also make certain assumptions on the parametric dependence of a, /, and L In particular, 
we shall suppose that, for some finite (preferably small) integer Q % a may be expressed as 


Q 

a (w t v\p) = ^ a q {p) a q {w 1 v) > Vw,ueX, V/iGD, 

( 7=1 


(2) 


for some o q : V — ► ]R and a q : X x X -> ]R, q = 1, . . . , Q. This “separability,” or “affine,” assumption 
on the par ametei dependence is crucial to computational efficiency; however, certain relaxations are 
possible — see Section 5.3.3. For simplicity of exposition, we assume that / and £ do not depend on p\ 
in actual practice, affine dependence is readily admitted. 

Our abstract problem statement is then: for any p £ V, find s(p) e R given by 


where u(p) € X is the solution of 


s{p) = £(u(p)) t 


( 3 ) 


a(it(/i), i;; jti) — /(u), V r £ X. (4) 

In the language of the introduction, a is our partial differential equation (in weak form), p is our 
parameter, u(p) is our field variable, and s(p) is our output. For simplicity of exposition, we may on 
occasion suppress the explicit dependence on p. 
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2.2 Particular Instantiations 

We indicate here a few instantiations of the abstract formulation; these will serve to illustrate the 
methods (for coercive, symmetric problems) of Sections 3 and 4. 


2.2.1 A Thermal Fin 



Figure 1: Two- and Three-Dimensional Thermal Fins. 


In this example we consider the two- and three-dimensional thermal fins shown in Figure 1; these 
examples may be (interactively) accessed on our web site 1 . The fins consist of a vertical central “post” 
of conductivity k 0 and four horizontal “subfins” of conductivity k\ i = The fins conduct 

heat from a prescribed uniform flux source q ,f at the root r roo t through the post and large-surface- 
area subfins to the surrounding flowing air; the latter is characterized by a sink temperature u 0 and 
prescribed heat transfer coefficient h. The physical model is simple conduction: the ternperatuie field 
in the fin, u, satisfies 

Y f HiVH f h{u-u 0 )v= [ q"v, (5) 

Yo Ah Ja n\f ( „ot 

where ft* is that part of the domain with conductivity k\ and dQ. denotes the boundary of fh 

We now (t) nondimensionalize the weak equations (5), and (n) apply a continuous piecewise- affine 
transformation from ft to a fixed (ju-independent) reference domain ft [15]. The abstract problem 
statement (4) is then recovered for fx = {L\ /e 2 , L 3 , L 4 , Bi, L, £}, V — [0.1,10.0] x [0.01,1.0] x 
[2.0, 3.0] x [0.1 x 0.5], and P ~ 7; here are the thermal conductivities of the “subfins” (see 

Figure 1) relative to the thermal conductivity of the fin base; Bi is a nondimensional form of the heat 
transfer coefficient; and, L, t are the length and thickness of each of the “subfins” relative to the length 
of the fin root i\oot- ft is readily verified that a is continuous, coercive, and symmetric; and that the 
“affine” assumption (2) obtains for Q ~ 16 (two-dimensional case) and Q - 25 (three-dimensional case). 
Note that the geometric variations are reflected, via the mapping, in the er^/i). 

For our output of interest, s(^), we consider the average temperature of the root of the fin nondi- 
mensionalized relative to q" , fc°, and the length of the fin root. This output may be expressed as 

1 Fin 2D: http://augustine.mit.edu/fin2d/fin2d.pdf and Fin 3D: http://augustine.mit.edu/fin3d-l/fin3d_l.pdf 
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HZ') - £("(/'))> "'here f(t>) - J r ^ v. It is readily shown that this output functional is bounded and 
also “compliant”: £(v) = f{v), Vu e X . 

2.2.2 A Truss Structure 



Figure 2: A Truss Structure 


We consider a prismatic microtruss structure [8, 27] shown in Figure 2; this example may be (interac- 
tively) accessed on our web site 2 . The truss consists of a frame (upper and lower faces, in dark gray) 
and a coie (tiusses and middle sheet, in light gray). The structure transmits a force per unit depth 
F u inform ly distributed over the tip of the middle sheet f 3 through the truss system to the fixed left 
wall IV The physical model is simple plane -strain (two-dimensional) linear elasticity: the displacement 
field Ui, i — 1,2, satisfies 


f On p 0u k 

In ' 3kl dx t 


wlieie Q is tlie truss domain, E tJ ki is the elasticity tensor, and X refers to the set of functions in H 1 (il) 
which vanish on Fq. We assume summation over repeated indices. 

We now (i) liondimensionalize the weak equations (6), and («) apply a continuous piccewise-affinc 
transformation from O to a fixed (/^-independent) reference domain ft. The abstract problem state- 
ment (4) is then recovered for /x = {t fi t u H, ^},D= [0.08, 1.0] x [0.2, 2.0] x [4.0, 10.0] x [30.0°, 60.0°], 
and P = 4; here t/ and t t are the thicknesses of the frame and trusses (normalized relative to t c ) y 
respectively; H is the total height of the microtruss (normalized relative to t c ); and 0 is the angle be- 
tween the tiusses and the faces. The Poisson’s ratio, v = 0.3, and the frame and core Young’s moduli, 
E f - 75 GPa and E c = 200 GPa, respectively, are held fixed. It is readily verified that a is continuous^ 
coercive, and symmetric; and that the “affine” assumption (2) obtains for Q = 44. 

Our outputs of interest are (z) the average downward deflection (compliance) at the core tip, r 3 , 
nondimensionalized by F/Ej\ and ( ii ) the average normal stress across the critical (yield) section 
denoted Tf in Figure 2. These compliance and noncoin pliance outputs can be expressed as s l {fx) ~ 
and s 2 (/i) = ^ 2 (u(/i)), respectively, where t l {v) = - v 2) and 



d X i p du k 
d Xj ijkl dxi 


aie bounded lineal’ functionals; here ** is any suitably smooth function in H l (Q s ) such that X ifii = 1 
on Tj and = 0 on where h is the unit normal. Note that s 1 (/x) is a compliant output, whereas 
s 2 (fx) is “noncompliant.” 

“ TRUSS: http : / / august ine . mit . edu/ simple jtruss/simple_truss . pdf 
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3 Reduced-Basis Approach 

We recall that in this section, as well as in Section 4, we assume that a is continuous, coercive, symmetric, 
and affine in (jl — see (2); and that ^(u) = /(u), which we denote “compliance.” 

3.1 Reduced-Basis Approximation 

We first introduce a sample in parameter space, Sn = {/>i> ■ • • , / 4 n} ; where m G V> i = 1, . . ■ , A; see 
Section 3.2.2 for a brief discussion of point distribution. We then define our Lagrangian [22] reduced- 
basis approximation space as IV/v = span (£ n = u(fi n ), n = 1, . . . , N] , where u(/A n ) € X is tlie solution 
to (4) for f.L = fi n . In actual practice, u(/* Tl ) is replaced by an appropriate finite element approximation 
on a suitably fine truth mesh; we shall discuss the associated computational implications in Section 3.3. 
Our reduced -basis approximation is then: for any {jl € T> ) find syv(/i) = £(uw(f.i)), where G 

is the solution of 

a(u^(/i), v; fi) — ^(u), V t/ e Wjv- (7) 

Non-Galcrkin projections are briefly described in Section 5.3.1. 


3.2 A Priori Convergence Theory 

3.2.1 Optimality 

We consider here the convergence rate of w(aO an< ^ s n(aO ^ 5 (m) as N — > oo. To begin, it is 

standard to demonstrate optimality of u^(^i) in the sense that 

inf ||u(/.i) - M/v||x ■ ( 8 ) 

w/v€ VV/sf 

(We note that, in the coercive case, stability of our (“conforming” ) discrete approximation is not an issue; 
the noncoercive case is decidedly more delicate (see Section 5.3.1).) Furthermore, for our compliance 
output, 

s(fj) = + £(u — un) = s/v(m) +o(u , u - un\p) — sw{ii) + a(u — UN, u (9) 

from symmetry and Galcrkin orthogonality. It follows that <s (/z) - 3 N (fj) converges as the square of the 
error in the best approximation and, from coercivity, that s^(/i) is a lower bound for s(/i). 


||ti(/i) — ii.w(pt)|| 



3.2,2 Best Approximation 

It now remains to bound the dependence of the error in the best approximation as a function of N. At 
present, the theory is restricted to the case in which P = 1, V = [0,^ max ], and 

a(w, v\ /i) = ao(w, t>) 4- /iai(iLP, u), (10) 

where a 0 is continuous, coercive, and symmetric, and ai is continuous, positive semi-definite (a-i > 

0, Vie G A), and symmetric. This model problem (10) is rather broadly relevant, for example to variable 
orthotropic conductivity, variable rectilinear geometry, variable piecewise-constant conductivity, and 
variable Robin boundary conditions. 

We now suppose that the jz n , n = 1, . . . , A, are logarithmically distributed in the sense that 

In (A (x ri + l) = ^ In (A \x max 4- l) , n = 1, . . . , N, (11) 

where A is an upper bound for the maximum eigenvalue of a\ relative to a 0 . (Note A is perforce bounded 
thanks to our assumption of continuity and coercivity; the possibility of a continuous spectrum does 
not, in practice, pose any problems.) We can then prove [18] that, for N > N cr j t = eln(A /i max T 1), 

inf r ||u(^) - ujn(a 4 )IIx <(14* Mmax A) ]|u(0)||x exp | ( jv crit _ ij | ’ V H € D ■ (12) 
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N 

|s(/i) - s A r(/i)|/s(/j) 

a a '(p)/s(p) 

V.win) 

10 

1.29 x lO 172 

8.60 x 10^ 2 

2.85 

20 

1.29 x 1CT 3 

9.36 x 10~ 3 

2.76 

30 

5.37 x 10~ 4 

4.25 x 10~ 3 

2.68 

40 

8.00 x 10- 5 

5.30 x 10~ 4 

2.86 

50 

3.97 x lO -5 

2.97 x 10“ 4 

2.72 

60 

1.34 x 10- 5 

1.27 x 10“ 4 

2.54 

70 

8.10 x 10-° 

7.72 x 10~ 5 

2,53 

80 

2.56 x lO -6 

2.24 x 10“ 5 

2.59 


Table 1: Error, error bound (Method I), and effectivity as a function of AT, at a particular representative 
point \l € Z), for the two-dimensional thermal fin problem (compliant output). 


We observe exponential convergence, uniformly (globally) for all p in with only very weak (loga- 
rithmic) dependence on the range of the parameter (/x lllax ). (Note the constants in (12) are for the 
particular case in which (•, )x ~ ao(-, •)*) 

The proof exploits a parameter-space (non- polynomial) interpolant as a surrogate for the Galerkin 
approximation. As a result, the bound is not always “sharp”: we observe many cases in which the 
Galerkin projection is considerably better than the associated interpolant; optimality (8) chooses to 
illuminate only ceitain points fi n , automatically selecting a best “sub— approximation” amongst all 
(combinatoiially many) possibilities we thus see why reduced— basis state-space approximation of 
s(/i) via u(/x) is preferred to simple parameter- space interpolation of s(p) (“connecting the dots”) via 
(p ni s(p n )) pairs. We note, however, that the logarithmic point distribution (11) implicated by our 
interpolant- based arguments is not simply an artifact of the proof: in numerous numerical tests, the 
logarithmic distribution performs considerably (and in many cases, provably) better than other more 
obvious candidates, in paiticular for large ranges of the parameter. Fortunately, the convergence rate 
is not too sensitive to point selection: the theory only requires a log “on the average” distribution [18]; 
and, in practice, A need not be a sharp upper bound. 

The result (12) is certainly tied to the particular form (10) and associated regularity of u(p). How- 
ever, we do observe similar exponential behavior for more general operators; and, most importantly, the 
exponential convergence rate degrades only very slowly with increasing parameter dimension, P. We 
present in Table 1 the error |sQx) - s N (fj,)\/s(fj.) as a function of TV, at a particular representative point 
/i in 72, foi the two-dimensional thermal fin problem of Section 2.2.1; we present similar data in Ta- 
ble 2 foi the tiuss problem of Section 2.2.2. In both cases, siiice tensor-product grids are prohibitively 
profligate as P incieasos, the /x n are chosen “log-randomly” over Zb we sample from a multivariate 
uniform probability density on log(^). We observe that, for both the thermal fin (P = 7) and truss 
(P = 4) problems, the error is remarkably small even for very small TV; and that, in both cases, very 
rapid convergence obtains as TV -> oo. W r c do not yet have any theory for P > 1. But certainly the 
Galerkin optimality plays a central role, automatically selecting “appropriate” scattered- data subsets 
of S/v and associated “good” weights so as to mitigate the curse of dimensionality as P increases; and 
the log-random point distribution is also important — for example, for the truss problem of Table 2, a 
non-logarithmic uniform random point distribution for S N yields errors which are larger by factors' of 
20 and 10 for TV = 30 and 80, respectively. 

3.3 Computational Procedure 

The theoretical and empirical results of Sections 3.1 and 3.2 suggest that TV may, indeed, be chosen 
very small. We now develop off-line/on-line computational procedures that exploit this dimension 
reduction. 
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N 

l«(/*) - sw(/0IMm) 

Ajv(aO/s(/j) 

'In (/>■) 

10 

3.26 x 10' 

6.47 x 10“'^ 

1.98 

20 

2.56 x 10 4 

4.74 x 10 4 

1.85 

30 

7.31 x 10~ 5 

1.38 x 10^ 4 

1.89 

40 

1.91 x 10~ 5 

3.59 x 10“ 5 

1.88 

50 

1.09 x 10~ 5 

2.08 X 10~ 5 

1.90 

60 

4.10 x 10~ 6 

8.19 x 10- G 

2.00 

70 

2.61 x 10~ 6 

5.22 x 10“ 6 

2.00 

80 

1.19 x 10~ 6 

2.39 x 10~ 6 

2.00 


Table 2: Error, error bound (Method II), and effectivity as a function of N, at a particular representative 
point p G V , for the truss problem (compliant output). 


We first express un(p) as 

N 

un{p) = X! UN iO- 4 ) 0 “ (^n(/0) 7 C» ( 13 ) 

i-l 

where u N (p) 6 R N ; we then choose for test functions v = i = 1, . . . , iV- Inserting these representa- 
tions into (7) yields the desired algebraic equations for u N (p) € R N , 

An(p) ^n(p) = £n- ( 14 ) 

in terms of which the output can then be evaluated as s N (p) = F^ u N (p). Here A N (p) G H NxN is 
the SPD matrix with entries An ij(fTj = u{CjXi\p)i 1 ^ h 3 Si N, and * s the ^ oac ^ ( a nd 

“output”) vector with entries Fn i = / (Ct ) 5 i = H • ■ • i N. 

We now invoke (2) to write 

Q 

A n ij{fx) = a(0, C*; M) = X] * V) C<) . ( 15 ) 

q-1 


AwM — Xv N ’ 

9" 1 

where the ^ € R Kx " are given by = a*(Cj,C0. * < hi < N > 1 S Q < Q- The off-line/on-line 

decomposition is now clear. In the off-line stage, we compute the u{p n ) and form the A? N and F N : this 
requires JV (expensive) “a” finite element solutions and 0(QN 2 ) finite-element- vector inner products. 
In the on-line stage, for any given new p, we first form A N from (15), then solve (14) for u N {y), and 
finally evaluate s N (p) = F% u N (fJi): this requires 0{QN 2 ) + 0(\N 2 ) operations and 0{QN 2 ) storage. 

Thus, as required, the incremental, or marginal, cost to evaluate s^(/r) for any given new p as 
proposed in a design, optimization, or inverse-problem context — is very small: first, because N is very 
small, typically 0(10) — thanks to the good convergence properties of Wn\ and second, because (14) 
can be very rapidly assembled and inverted — thanks to the off-line/on-line decomposition (see [5] for 
an earlier application of this strategy within the reduced-basis context). For the problems discussed 
in this paper, the resulting computational savings relative to standard (well-designed) finite-element 
approaches are significant — at least 0(10), typically 0(100), and often 0(1000) or more. 

4 A Posteriori Error Estimation: Output Bounds 

From Section 3 we know that, in theory, we can obtain s N {p) very inexpensively: the on-line compu- 
tational effort scales as 0(|iV 3 ) + 0{QN 2 )\ and N can, in theory , be chosen quite small. However, in 
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practice, we do not know how small N can be chosen: this will depend on the desired accuracy, the 
selected output(s) of interest, and the particular problem in question; in some cases N = 5 may suffice 
while in other cases, N = 100 may still be insufficient. In the face of this uncertainty, either too many 
or too few basis functions will be retained: the former results in computational inefficiency; the latter 
in unacceptable uncertainty particularly egregious in the decision contexts in which reduced- basis 
methods typically serve. We thus need a posteriori error estimators for s N . Surprisingly, a posteri- 
ori error estimation has received relatively little attention within the reduced basis framework [19] 
even though reduced-basis methods are particularly in need of accuracy assessment: the spaces are 
ad hoc and pre- asymptotic, thus admitting relatively little intuition, “rules of thumb,” or standard 
approximation notions. 

Recall that, in this section, we continue to assume that a is coercive and symmetric, and that i is 
“compliant.” 

4.1 Method I 

The approach described in this section is a particular instance of a general “variational” framework for 
a posteriori error estimation of outputs of interest. However, the reduced-basis instantiation described 
here differs significantly from earlier applications to finite element discretization error [16, 14] and 
iterative solution error [20] both in the choice of (energy) relaxation and in the associated computational 
artifice. 


4.1.1 Formulation 

We assume that we are given a positive function y(p) \ V — > ]R +1 and a continuous, coercive, symmetric 
(/x-independent) bilinear form a : X x X -> IR, such that 

^olMIx < < a(v^v\p)^ V u 6 X, V /x e D (16) 

for some positive real constant a 0 . We then find e(ji) G X such that 

9 {^) a(e(p) } v) = R(v;u N (p)‘f.i) ) Vu G X, (17) 

where for a given w G X, R(v\w\fi) = C(v) - a(w,v\n) is the weak form of the residual. Our lower and 
upper output estimators are then evaluated as 

s N {p) = s N (p), and s%{p) = s N {p) + (18) 

respectively, where 

S g{p) a(e(n),e{n)) (19) 

is the estimator gap. 

4.1.2 Properties 

We shall prove in this section that s~ N {p) < s{p) < s%(p), and hence that | s (/i) - s N ( p )| = s (p) - 
s N (p.) < A N (p). Our lower and upper output estimators are thus lower and upper output bounds; and 
our output estimator gap is thus an output bound gap - — a rigorous bound for the error in the output 
ol interest. It is also critical that A/v(p) be a relatively sharp bound for the true error: a poor (overly 
laige) bound will encourage us to refine an approximation which is, in fact, already adequate — with 
a corresponding (unnecessary) increase in off-line and on-line computational effort. We shall prove in 
this section that A N (p.) < ^(s(p) - s N (p)), where 70 and Oq are the N -independent a-continuity and 
<?(/t)a-coercivity constants defined earlier. Our two results of this section can thus be summarized as 


1 < Vn(p) < Const, VAT, 


( 20 ) 
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where 


VN O) 


A,v(m) 

- 5 n(/0 


(21) 


is the efjectivity , and Const is a constant independent of iV. We shall denote the left (bounding 
property) and right (sharpness property) inequalities of (20) as the lower effcctivity and upper effectivity 
inequalities, respectively. 

We first prove the lower effectivity inequality (bounding property): s^(p) < s(fi) < (ft) , V/i € X>, 
for sjf(fi) and s£(/a) defined in (18). The lower bound property follows directly from Section 3.2.1. To 
prove the upper bound property, we first observe that R(v\ ujv; fj.) ~ a(u({L)-u N (ii) } v] fi) = a(e(p), /i), 

where e(^) = u(/*) - un(p)\ we may thus rewrite (17) as f/(/i)a(e(/i), v) = a(e(/i), u; /a) V u 6 X. We 
thus obtain 


4(/i)ft(e, e) = <?(//)a(e - e, e - e) + 2 g{fi)a(e ) e) - g{^)a{e } e) 

= 9(/j)a(e - e, e - e) -f (a(e, e; //) - <K/^)a(e, e)) + a(e, e; //) 

> a(e,e;/A) (22) 


since y(^) a(e(/x) - e(/^),e(/i.) - e(/i)) > 0 and a(e(/i), e(^); p) - p(/i) d(e(//), c(/a)) > 0 from (16). 
Invoking (9) and (22), we then obtain s(/a) - s N (f.t) — a(e(/i), c(/a);/a) < p(/i) a(e(fi),e(fi))\ and thus 
5(/i) < sjv(/i) + 0M e(/x)) = s£(/z), as desired. 

We next prove the upper effectivity inequality (sharpness property): 


Vn(p) = 


_ < 2°, v tv. 

s(^)-sn(m) «o’ 


To begin, we appeal to a-continuity and </(/;)d-coercivity to obtain 


a(e{ij,),e{n);ii) < 


7o 3 (m) - 

a 

«0 


(e(/i),e(/x)). 


(23) 


But from the modified error equation (17) we know that g(/r)d(e(jt), e(p)) = R{e(f l);h) = a(e(/x), e(^i); ^.). 
Invoking the Cauchy-Schwartz inequality, we obtain 


g{n)a(e,e) 


a(e,e;/i) < («(e,e;/i)) l/2 (a(e,e;^t)) 1/2 < 



(g{fi) a(e,e)) 1/2 (a(e,e;/t)) 1/2 ; 


the desired result then directly follows from (19) and (9). 

We now provide empirical evidence for (20). In particular, we present in Table 1 the bound gap and 
effect! vities for the thermal fin example. Clearly is always greater than unity for any TV, and 

bounded — indeed, quite close to unity — as N -* oo. 


4.1.3 Computational Procedure 

Finally, we turn to the computational artifice by which we can efficiently compute A at (/a) in the on -line 
stage of our procedure. We again exploit the affine parameter dependence, but now in a less transparent 
fashion. To begin, we rewrite the “modified” error equation, (17), as 


a(e{fj,),v) 


553 |,w 


Q N 
<J=1 7 = 1 


v v e X, 


where we have appealed to our reduced-basis approximation (13) and the affine decomposition (2). It 
is immediately clear from linear superposition that we can express e(fi) as 




{h)un j(n)Z] 


(24) 
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H) 


where i 0 £ X satisfies a(zo,v) = C(v), V v £ X, and Z] £ X,j = 1 , TV, q = 1 Q, satisfies 

a{Zj,v) = -a^Cj.u), V ti £ A'. Inserting (24) into our expression for the upper bound, s%(n) 
sn(p) + ^(/i)a(fi(/z), e(/i)), we obtain 


e + 


>v(p) = + 




Q N Q Q N N 

C ° +2 X^Xl a<7 (^)^ v >(^)A j + X] TV* (/*)<** (^)'^iVi(M) u N: 

<j--l y' = 1 j = l j' = l 


<?=1 j = l 


■0*>r& 


where c 0 = a(i 0 i ^o), Aj = a(*o, £?), and TJV' = a(zj } £?'). 


(25) 


The off-line/ on line decomposition should now be clear. In the off-line stage we compute zq and 
z], j = 1, . . . , N, q = 1, . . . , Q | and then form c 0 , AJ, and this requires QA/’ + 1 (expensive) “a” 
finite element solutions, and O(Q^iV^) finite-element-vector inner products. In the on-line stage, for 
any given new we evaluate as expressed in (25): this requires 0(Q 2 N 2 ) operations and 0(Q 2 N 2 ) 

storage (for c 0 , A*-, and TJJ,). As for the computation of sjv(/x), the marginal cost for the computation 

of s^(/i) for any given new (.1 is quite small — in particular, it is independent of the dimension of the 
truth finite element approximation space X. 

There are a variety of ways in which the off-line/on- -line decomposition and output error bounds 
can be exploited. A particularly attractive mode incorporates the error bounds into an on line adaptive 
process, in which we successively approximate s N (ff) on a sequence of approximation spaces W N >_ c 


Wn , N j = Nu2* — for example, may contain the N l - sample points of Syy closest to the new /i of 
interest - - until A;v' is less than a specified error tolerance. This procedure both minimizes the on-line 
computational effort and reduces conditioning problems — while simultaneously ensuring accuracy and 
certainty. 

The essential advantage of the approach described in this section is the guarantee of rigorous bounds. 


There are, however, certain disadvantages. The first set of disadvantages relates to the choice of g(p) and 
a. In many cases, simple inspection suffices: for example, in our thermal fin problem of Section 2.2.1, 
g(f.t) = rnin^^i . q ^ 9 (/i) and a(w t v) = a q (w,v) yields the very good effectivities summarized 

in Table 1. In other cases, however, there is no self-evident (or readily computed [17]) good choice: for 
example, for the truss problem of Section 2.2.2, the existence of almost-pure rotations renders g(p) very 
small relative to y(p), with corresponding detriment to g The second set of disadvantages relates 
to the computational expense — the 0(Q) off-line and the 0(Q 2 ) on-line scaling induced by (24) and 
(25), respectively. Both of these disadvantages are eliminated in the “Method IF’ to be discussed in 
the next section; however “Method II” only provides asymptotic bounds as N — > 00 . The choice thus 
depends on the relative importance of absolute certainty and computational efficiency. 


4.2 Method II 

As already indicated, Method I has certain limitations; we discuss here a Method II which addresses 
these limitations — - albeit at the loss of complete certainty. 


4.2.1 Formulation 

To begin, we set M > N } and introduce a parameter sample Sm = {Mil ♦ * • , y>M } and associated 
reduced-basis approximation space W M = span {< m = ti(/i m ), m = 1, . . . , M} ; for both theoretical and 
practical reasons we require S N C Sm and therefore W N C W M - The procedure is very simple: we first 
find UA/(/i) 6 Wm such that a(ujv,/(/i), v\ p) = /(u),V v e Wm\ we then evaluate sa/(aO = ^(um(/x)); 
and, finally, we compute our upper and lower output estimators as 

s n,m(p) = s n,mM = s n(v) + (26) 

where the estimator bound gap, is given by 

&N,\i(n) = I (sm(p) - SN{fl)) 

T 


( 27 ) 
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for some r G (0, 1). The effectivity of the approximation is defined as 


>7n,a/(a0 “ 


A/y^/pi) 

s(/i) - S N {p) 


(28) 


For our purposes here, we shall consider M — 2N . 


4.2.2 Properties 

As for Method I, we would like to prove the effectivity inequality 1 < ?7/v,2/v(/0 < Const, ViV. However, 
we will only be able to demonstrate an asymptotic form of this inequality. Furthermore, the latter shall 
require — and we shall make — the hypothesis that 




as N -* oo 

s(fl) - S N (fl) 


(29) 


We note that the assumption (29) is certainly plausible: if our a priori bound of (12) in fact reflects 
asymptotic behavior, then s(/i) — sjsj(fi) ~ c\e~ C2 ^ , s(p) — S2 n{t) ~ c i e 2c2N y and hence £/v,2n(m) ~ 
e“ CaAr , as desired. 

We first prove the lower effectivity inequality (bounding property): s^ 2N (fx) < s(fi) < s^ 2N (fi) y as 
N —> oo. To demonstrate the lower bound we again appeal to (9) and the cocrcivity of a; indeed, this 
result (still) obtains for all N. To demonstrate the upper bound, we write 


5 n,27v(m) 


s + (d - l) (s(m) - sjv(m)) - To 1 ) - s 2 n{i*)) 

8 + (”[1 ~~ £ N ,2N (m)] “ l) ( 5 (f0 — s N(l ^))- 


(30) 

(31) 


We now recall that s(fi) - s,y(m) > 0, and that 0 < r < 1 — that is, 1/r > 1; it then follows from (31) 
and our hypothesis (29) that there exists a finite N such that 

4,2 *(/i)-*(M)>0, V7V> N\ (32) 


This concludes the proof: we obtain asymptotic bounds. 

We now prove the upper effectivity inequality (sharpness property) . From the definitions of rfw^N (m)> 

A/v, 2 n(a 0 and £ n,2n(m)> we directly obtain 


1 S2n{^) ~ SN(M) 

1 («2n(/0 “ 5 (A)) ~ ( 5 7V(M) “ s (p)) 

(33) 

r s(n) - s N (fJ.) 

T - S N (ll)) 

= -(1 ~ £N,2n{v))- 

t" 


(34) 


It is readily shown that ;;n,2n(m) bounded from above by 1/r for all N : we know from (9) that 
e N2 N{f J ') is strictly non-negative. It can also readily be shown that r}s t 2N(p) is non-negative: since 
c it follows from (8) (for (•, -)x = a(', *; a 0) and (9) that s(/r) > S2 n(^) > $n(m)> an d hence 
cn, 2 n(h) < 1- We thus conclude that 0 < w, 2 n(m) < V r for a11 N - Furthermore, from our hypothesis 
on £/v,2n(m)> (29), we know that VN,2N{p) will tend to 1/r as N increases. 

The essential approximation enabler is exponential convergence: we obtain bounds even for rather 
small N and relatively large r. We thus achieve both “near” certainty and good effect ivi ties. We 
demonstrate this claim in Table 2, in which we present the bound gap and effectivity for our truss 
example of Section 2.2.2; the results tabulated correspond to the choice r = 1/2. We clearly obtaift 
bounds for all N\ and we observe that 77 /v,2n(a) does, indeed, rather quickly approach 1/r. 
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4.2.3 Computational Procedure 

Since the error bounds are based entirely on evaluation of the output, we can directly adapt the off 
line/on-line procedure of Section 3.3, Note that the calculation of the output approximation .s \r(/i) 
and the output bounds are now integrated: A N (fi) and F N (p) (yielding s jV (V)) are a sub-matrix and 
sub-vector of A 2JV (/i) and F 2N (}l) (yielding a* 2 n(/x), A^/o /vQu), and s^- w (/i)), respectively. In the 
off-line stage, we compute the u(^ Tl ) and form the and £ 2N : this requires 2N (expensive) V’ 
bnite element solutions, and 0(4QN 2 ) finite-element-vector inner products. In the on-line stage, for 
any given new /i, we first form A N {p) t F N and A 2N (ff), F 2N , then solve for u N {ix) and u 2N {n) t and 
finally evaluate s% 2N {ff): this requires 0{AQN 2 ) + 0(~N Z ) operations and 0(4QIV 2 ) storage. The 
on- line effort for this Method II predictor/error estimator procedure (based on and s 2 /v(/a)) will 

thus require eightfold more operations than the “predictor-only” procedure of Section 3. 

Method II is in some sense very naive: we simply replace the true output s(fi) with a filler- 
approximation surrogate s 2 n(m)- (There are more obscure ways to describe the method — in terms of a 
reduced -basis approximation for the error — however there is little to be gained from these alternative 
iuteipietations.) ffhe essential computation enabler is again exponential convergence, which permits 
us to choose M = 2 N — hence controlling the additional computational effort attributable to error 
estimation - - while simultaneously ensuring that e Ny2 N{ff) tends rapidly to zero. Exponential conver- 
gence also ensures that the cost to compute both s N (p) and s 2 /v(/j) is “negligible.” In actual practice, 
since s 2 w(^) is available, we can of course take s 2 /v(m), rather than s/v(/i), as our output prediction; 
this gieatly impioves not only accuracy, but also certainty - — A/v i2 /v(m) is almost surely a bound for 
- S 2 /v(m), albeit an exponentially conservative bound as N tends to infinity. 

5 Extensions 

5.1 Noncompliant Outputs and Nonsymmetric Operators 

In Sections 3 and 4 we formulate the reduced- basis method and associated error estimation procedure 
for the case of compliant outputs, i(v) = /( t>), Vu 6 X . We briefly summarize here the formulation 
and theory for more general linear bounded output functionals; moreover, the assumption of symmetry 
(but not yet coercivity) is relaxed, permitting treatment of a wider class of problems — a representative 
example is the convection-diffusion equation, in which the presence of the convective term renders the 
operator nonsymmetric. We first present the reduced-basis approximation, now involving a dual or 
adjoint problem; we then formulate the associated a posteriori error estimators; and we conclude with 
a few illustrative results. 

As a preliminary, we first generalize the abstract formulation of Section 2.1. As before, we define the 
“primal” problem as in (4), however we of course no longer require symmetry. But we also introduce 
an associated adjoint or “dual” problem: for any X, find 6 X such that 

a(v,i Viiel; (35) 

recall that £{v) is our output functional. 

5.1.1 Reduced-Basis Approximation 

To develop the reduced-basis space, we first choose — randomly or log-randomly as described in Sec- 
tion 3.2 — a sample set in parameter space, S N/2 = {/z x , . . . , /x*r/ 2 }, where m £ V } i = l^.^N/2 
( N even); we next define an “integrated” Lagrangian reduced-basis approximation space, = 
span {(u(/x n ),^(/i„)), 1 ,..., N/2}. 

For any fj £ D, our reduced basis approximation is then obtained by standard Galerkin projection 
onto W N (though for highly nonsymmetric operators minimum residual and Petrov- Galerkin projections 
are attractive — stabler — alternatives). To wit, for the primal problem, we find u N (fi) e W N such 
that a(upj (/i), v] p) = f(v ), Vv £ W/v; and for the adjoint problem, we define (though, unless otherwise 
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indicated, do not compute) Vn{i i ) £ Wn such that a(t>, p) — ^( v )> ^ The ice! need- 

basis output approximation is then calculated from s^(p) = £{un(p))- 

Turning now to the a priori theory, it follows from standard arguments that un{p) and ip N {p) arc 
“optimal” in the sense that 

||w(M) - U N (jl)\\x < (l+^) 

- *n(h)u < ( i + U1 w » i w 

The best approximation analysis is then similar to that presented in Section 3.2. As regards our output, 
we now have 

\s{p)-s N (p)\ = |£(u(/i))-^(ujv(p))| - \a{u-u Ni fa v)\ = \a(u-u Ni tM/w p)\ < 'yb||*u-n w ||x||^-^Af|U 

(3G) 

from Galerkin orthogonality, the definition of the primal and the adjoint problems, and the Cauchy- 
Schwartz inequality. We now understand why we include the ipipn) in Wjsfi to ensure that ||?/>(/x) — 
V’nOOIU is small. We thus recover the “square” effect in the convergence rate of the output, albeit 
(and unlike the symmetric case) at the expense of some additional computational effort — the inclusion 
of the ip(}in) in Wn\ typically, even for the very rapidly convergent reduced-basis approximation, the 
“fixed error-minimum cost” criterion favors the adjoint enrichment. 

For simplicity of exposition (and to a certain extent, implementation), we present here the “inte- 
grated” primal-dual approximation space. However, there are significant computational and condition- 
ing advantages associated with a “non-integrated” approach, in which we introduce separate primal 
(u(p n )) and dual (i>{p n )) approximation spaces for u(p) and i>(p), respectively. Note in the “non- 
integrated” case we are obliged to compute ^/v(/x)> since to preserve the output error “square effect” 
we must modify our predictor with a residual correction, f(ip^(p)) - a(u^(p), ^Pn(p)] p) (17]. Both the 
“integrated” and “non-integrated” approaches admit an off-line/ on-line decomposition similar to that 
described in Section 3.3 for the compliant, symmetric problem; as before, the on-line complexity and 
storage are independent of the dimension of the very fine (“truth”) finite element approximation. 


5.1.2 Method I A Posteriori Error Estimators 

We extend here the method developed in Section 4.1.2 to the more general case of noncompliant and 
nonsymmetric problems. We begin with the formulation. 

We first find e pr (/x) 6 X such that 

g(p) a(e pT (p),v) = jR pr (t>; tx/v(/x); /x), V v £ X, 

where R pr {v\w\p) = f(v) - a(tu,u;/x), V v £ X; and e du (/x) G X such that 

g(p)a(e du (p),v) = R 6u (v;^n{p)\ /x), V v G X, 


where R d "(v;w;p) = -t[y) - a(v,w\p),\f v G X. We then define 

**(/») = and (37) 

a n (m) = ^[a(« pr (M).s pr (M))]*[ a (« d "(M).« du 0*))] i - ( 38 ) 


Finally, we evaluate our lower and upper estimators as s^(p) = sn(/ 0 ± A w (/x). Note that, as before, 
g(fi) and a still satisfy (16); and that, furthermore, (16) will only involve the symmetric part of a. We 
define the effectivity as 

, X _ A w (a*) 

\s{n) - sn{i$\' 


(39) 
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note that s(}.i) — sjv(f.i) now has no definite sign. 

We now piove that our error estimators are bounds (the lower effectivity inequality): sjj(fi) < 

sb‘) < 40‘), V Af. To begin, we define e ± ( M ) = e pi (,r) rp ±e du (/r), and note that, from the cocrcivity 
of d, 


«5(/^)a(e pr - ^e ± ,e pr - ie*) = H.g(n)a(e pr , e pr ) + - ^(/i)a(e ± ,e p ') > 0, (40) 

where e p, (p) = u(fi) - u w (/i), e du (M) = 4>{p) - and n is a positive real number. Rom tlie 

definition of e^/r) and e pr (p), e du (p), we can express the “cross-term” as 


y(A i )“(e ± , e pr ) = /i p ‘(e pr ;u/v;^) T -R Au {e vr \^ N \tA = a(e pr ,e pr ;;u) qp -a(e pr ,e dl, ;p) 


a(e pl ,e pr ;p) ± -(s(p) - s/v(/r)), 


(41) 


since ii p, (e pr ;u N ;/i) = a(u,e p, » - a(u w ,e pr ;/r) = a(e pr , e pr ; ^), ii du (e pr ; V> w ; /t) = a(e pr ,^;^) - 
“( el if /hv,p) = a(e l ,e du ;p), and £(/*) — 1 ( un ) = — a(u — = — a(u — (by 

Galerkin orthogonality) = -o(e pr ,e d,, ;/i). We then substitute (41) into (40) to obtain 

±(s(p) - s lV (^)) < — fc(a(e pr ,e pi ;/i) - g(ji)a(e pr ,e pr )) + —^(e*,^) < ^^(e*,-!*), 

since k > 0 and a(e pr (/i), e pr (p); p) - fl (/i)a(e pr (p), e pr (/a)) > 0 from (16). 

Expanding e ± (/j) = e pr (^) =p £e du (/i) then gives 


±( s (^) - sn(/j.)) < 


p(m) 


«a(e pr ,e pr ) + ia(e du , e du ) ^ 2a(e pr , e du ) 


Ut 

± (s(p) - (s w (p) - ^a(e pr ,e du ))j < ^a(e pr ,e' ,r ) 
We now choose n(f.i) as 

^d(e d “(^),e du (p))V 


+ ^a(e du ,e d “). 


(42) 


«(/i) 


d(e pr (/i),e pr (p)) / 


so as to minimize the right-hand side (42); we then obtain 


MaO - s w (/t)| < A w (/i), (43) 

and hence sjO(p) 5- s (m) < s ti(^)- 

We now turn to the upper effectivity inequality (sharpness property). If the primal and dual errors 
are a-orthogonal, or become increasingly orthogonal as N increases, then the effectivity will not, in 
fact, be bounded as IV — » oo. However, if we make the (plausible) hypothesis that I s I n ) — sn(u)\ > 
C||e pr (Ai)||x||e du (/i)||x, then it is simple to demonstrate that 


win) < 


7o 

2 C Og 


(44) 


In particular, it is an easy matter to demonstrate that p 1/2 (p) (a(e pr (/i),e pr (/r))) 1/2 < -7"_|| e pr (^)|| x 

(note we lose a factor of 7 0 / relative to the symmetric case); similarly, g l/2 (fi) (a (e du (/i), e du (/i ))) 1/2 < 
11^ (/-Ollx* The desired lesult then directly follows from the definition of A^v(/i) and our hypothesis 
on \s(fi) - s n {im)\. 

Finally, tinning to computational issues, we note that the off-line/on-line decomposition described 
in Section 4.1 foi compliant symmetric problems directly extends to the noncompliant, nonsymrnetric 
case — except that we must compute the norm of both the primal and dual “modified errors,” with a 
concomitant doubling of computational effort. 


5 EXTENSIONS 


15 


5.1.3 Method II A Posteriori Error Estimators 

Wc discuss here the extension of Method II of Section 4.2 to noncompliant outputs and nonsymmetric 
operators. 

To begin, we set M > N, M even, and introduce a parameter sample S M / 2 = {/M 1 • • * » ^Af/ 2 } arid as- 
sociated “integrated” reduced-basis approximation space Wm = span { u(/j m ), V' (Mm), rn = 1, . . . , A'//2}. 
We first find -u A /(m) 6 IV* f such that = f(v), Vu € W M \ we then evaluate s M (/t) = 

f(u.A/(/i))' al 'd finally, we compute our upper and lower output estimators as 

• s W,m(/‘) = s N(/r) + - Sn) ± 

An.wW = -|sm(m) - sn(p)Ii ( 46 ) 

7 


for r € (0, 1). The effectivity of the approximation is defined as 


W,a/(m) 


KaO ~ sn(//)I' 


(47) 


We shall again only consider M — 2 TV. 

As in Section 4.2, we would like to prove that 1 < r/iv^yv (/^) < Const for sufficiently large N ; and, as 
in Section 4.2, we must again make the hypothesis (29). We first consider the lower effectivity inequality 
(bounding property), and prove that 


5 n,2n(^) — 5 (a0 - s n,2n0 1 )> &£ N * oo. 


(48) 


In particular, simple algebraic manipulations yield 


S N,2N (t 1 ) 
s N,2n(^) 


s(n) - — — \sn(v) ~ $2N (m)I x 

-L “ £N,2N 

s(/i) + |sn(aO “ s 2n(aOI X 

1 ~ 


“(1 “ £N, 2 n) “ 1 
“(1 - £N,2N) ~ 1 


52 n(/*) > 5 AT (A) 

* 2 n(m) < ’ 

S 2 f/(/*) > 5N(A0 
52N(M) < S N(M)‘ 


(49) 

( 50 ) 


The desired result then directly follows from our hypothesis on £jv t 2 JV» (29), and the range of r. 

The proof for the upper effectivity inequality (sharpness property) parallels the derivation of Sec- 
tion 4.2.2. In particular, we write 


*?n, 2 n(m) 


d|s 2 AT — S N \ __ ^|s 2 7V — s + S — S/vl 
1 5 Syv | |s Sj\r | 

-|1 — £N,2N\\ 

T 


(51) 


(52) 


from our hypothesis (29) we may thus conclude that 77 /v,2n(m) — ► 7 as N — * 00 . Note in the noncoin- 
pliant, nonsymmetric case we can make no stronger statement. 

We demonstrate our effectivity claims in Table 3, in which we present the error, bound gap, and 
effectivity for the noncompliant output (s 2 (/x), average stress) of the truss example of Section 2.2.2; 
the results tabulated correspond to the choice r = 1 / 2 . We clearly obtain bounds for all N\ and the 
effectivity rather quickly approaches 1 jr (for N > 120 , remains fixed at 1 /r = 2 . 0 ). 


5.2 Eigenvalue Problems 

We next consider the extension of our approach to symmetric positive definite eigenvalue problems. 
The eigenvalues of appropriately defined partial-differential-equation eigenproblems convey critical in- 
formation about a physical system; in linear elasticity, the critical buckling load; in dynamic analysis 
of structures, the resonant modes; in conduction heat transfer, the equilibrium timescales. Solution 
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N 

K/-0 - s N (fx)\/s(n) 

X,V,2 N(ix)/ S (^t) 

VN,2n(h) 

20 

2.35 x 10" a 

4.67 x 10“ 2 

1.99 

40 

1.74 x 10 -4 

3.19 x 10’ 4 

1.83 

60 

5.59 x 10- 5 

1.06 x 10~ 4 

1.90 

80 

1.44 x 10 -5 

2.73 x 10" s 

1.89 

100 

7.45 x 1(T 6 

1.40 x 10- 5 

1.88 


Table 3: Error, error bound (Method II), and effectivity as a function of iV, at a particular representative 
point /t € for the truss problem (noncompliant output). 


of large-scale eigenvalue problems is computationally intensive: the reduced-basis method is thus very 
attractive. 

The abstract statement of our eigenvalue problem is: find {ui(fi), A,(/j)) £ X x IR, i — 1, . . such 
that 

a(ui(fi),v\/.L) = A Vu 6 X, and = 1. (53) 

Here a is the continuous, coercive, symmetric form introduced earlier, and m is (say) the L 2 inner 
product over fl. The assumptions on a and rn imply the eigenvalues A ;(/i) will be real and positive. 
We order the eigenvalues (and corresponding eigenfunctions u t ) such that 0 < A^/i) < A 2 (fj.) < . . . ; we 
shall assume that and A 2 (/i) are distinct. We suppose that our output of interest is the minimum 

eigenvalue, 

= Ai(ju); (54) 

other outputs may also be considered. 

Following [11), we present here a reduced-basis predictor and a Method I error estimator for sym- 
metric positive-definite eigenvalue problems; we also briefly describe the simpler Method II approach. 

5.2.1 Reduced-Basis Approximation 

Wc sample - - randomly or log-randomly — our design space V to create the parameter sample S N/2 = 
{/U, . . . we then introduce the reduced-basis space W N = span{ui(/ii), u 2 (ni ),..., (/*w/ 2 ), 

^(^Ar/a)}. where we recall that and are the eigenfunctions associated with the first (small- 

est) and second eigenvalues X x (n) and A 2 (^), respectively. Note that has good approximation 
properties both for the first and second lowest eigenfunctions, and hence eigenvalues; this is required 
by the Method I error estimator to be presented below. Our reduced-order approximation is then: find 
(um A Ni(fi)) € Wn x IR, i — 1, . . . , N, such that 

“(ujv«(M)»v;/i) = Vu e W N) and m{u N ^(/i), u N <(/x)) = 1; (55) 

the output approximation is then s//(/i) = ^/j). 

The formulation admits an on-line/off-line decomposition [11] very similar to the approach described 
for equilibrium problems in Section 3. 

5.2.2 Method I A Posteriori Error Estimators 

As before, we assume that we are given a positive function g(fi) : T) — > JR+ and a continuous, coercive, 
symmetric bilinear form a(iu,u) : X x X -» IR, that satisfy the inequality (16). We then find e(/u) € X 
such that 

= IV im(u w - a(u w v;/j)] , Vu 6 X, (56) 

in which the right-hand side is the eigenproblem equivalent of the residual. We then evaluate our 
estimators as 


— 1 (a^) » 


s wW = Ajy i(^) - A/v(/i), 
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N 

I-Mm) - W/‘)l/M/‘) 

A,v(m)/ Al(m) 

Wv(/*) 

10 

1.19 x 10 - " 2 

0.60 x 10- 2 ' 

5.63 

20 

i.08 x nr 3 

7.19 x 10~ 3 

6.65 

30 

0.20 x io~ 4 

3.19 x IO' 3 

5.17 

40 ; 

1.72 x 10~ 4 

1.55 x 10~ 3 

9.44 

50 

3.47 x 10~ 5 

4.00 x 10 -4 

11.74 


Table 4: Error, error bound (Method I), and effectivities as a function of N } at a particular representative 
point ft G D, for the thermal fin cigenproblen. 


Aj v(^) = -^~“a(e(/r),e(/t)), 

™(/x) 

where 5(jj.) = 1 - and r 6 (0, 1). The effectivity is defined as ,00 -a, (,.))• 

We now consider the lower and upper effectivity inequalities. As regards the lower effectivity in- 
equality (bounding property), we of course obtain s^(/i) > Ai^t), V N . The difficult result is the lower 
bound: it can be proven [11] that there exists an N*{S N / 2) h) such that s n (ja) < Ai(/.t), ViV > NT 
In practice, N* — due to the good (theoretically motivated) choice for <5(/i); there is thus very little 
uncertainty in our (asymptotic) bounds. We also prove in [11] a result related to the upper effectiv- 
ity inequality (sharpness property); in, practice, very good effectivities are obtained. To demonstrate 
these claims we consider the eigenvalue problem associated with (the homogenous version) of our two- 
dimensional thermal fin example of Section 2.2.1. We present in Table 5.2.2 the error, error bound, and 
effectivity as a function of JV at a particular point fi G V. We observe rapid convergence, bounds for 
all N } and good effectivities. 

Finally, we note that our output estimator admits an off-line/on-line decomposition similar to 
that for equilibrium problems; the additional terms in (56) are readily treated through our affine 
expansion/ linear superposition procedure. 


5.2.3 Method II A Posteriori Error Estimators 

For Method II, we no longer require an estimate for the second eigenvalue. We may thus define S u = 
{^i,. 1>3 jz n }, W n = span{ui(/ii), i = and (for M = 2 N) S 2N - {m» • • ■ >^2/v} 3 S Ny 

W 2N — span(‘Ui(/it),« = 1, . . . ,2A r } D W The reduced basis approximation now takes the form (53), 
yielding s N (n) = A n i(m) and (for N 2 N) s 2N (v) = A 2N i(/ 0- 0ur estimators are then given by 


s n,2n(^) = = A N i(/i) - 

Ayv,2 aKaO " ~( s n(aO - 5 2 n(^)) ( 57 ) 

r 


for r G (0,1). The effectivity 77n, 2A'(^) is defined as for Method I. 

For the lower effectivity inequality (bounding property), wc of course retain s^ 2N {(^) > A x (pt), VN. 
We also readily derive sJf t7N (fi) = Ai - (Aat i ~ Ai)(£(l -e N%2N ) - 1); under our hypothesis (29), we thus 
obtain asymptotic bounds as N -* oo. For the upper effectivity inequality (sharpness property), we 
directly obtain = £(1 -£n,2n)- By variational arguments it is readily shown that 0 < £/v,2 n < h 

we thus conclude that 0 < VN, 2 N ^ VAT. Additionally, under hypothesis (29), we deduce that 
Vn,2N ^ as A ^ oo. 


5.3 Further Generalizations 

In this section we briefly describe several additional extensions of the methodology. In each case we 
focus on the essential new ingredient; further details (in most cases) may be found in the referenced 
literature. 
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5.3.1 Noncoercive Linear Operators 

The archetypical noncoercive linear equation is the Helmholtz, or reduced-wave, equation; many (e.g., 
hivexse scattering) applications of this equation arise, for example, in acoustics and electromagnetics. 
The essential new mathematical ingredient is the loss of coercivity of a. In particular, well -posed ness 
is i low ensured only by the inf -sup condition: there exists positive P 0i /?(/*), such that 


0 < Po 




(58) 


Two numerical difficulties arise due to this “weaker” stability condition. 

Hie first difficulty is preservation of the inf-sup stability condition for finite dimensional approxima- 
tion spaces, do wit, although in the coercive case restriction to the space W p actually increases stability, 
in the noncocicive case lcstiictiou to the space VF/y can easily decrease stability: the relevant supremiz- 
ers may not be adequately represented. Loss of stability can, in turn, lead to poor approximations - the 
inf-sup parameter enters in the denominator of the a priori convergence result. The second numerical 
difficulty is estimation of the inf-sup parameter, which for noncoercive problems plays the role of g(p,) 
iu Method I a posteriori error estimation techniques. In particular, P(}i) can not typically be deduced 
analytically, and thus must be evaluated (via an eigenvalue formulation) as part of the reduced-basis 
approximation. Our resolution of both these difficulties involves two elements [17]: first, we consider 
projections other than standard Galerkin; and second, we consider “enriched” approximation spaces. 

In one approach [17], we pursue a minimum-residual projection: the (low-dimensional) infimizing 
space contains both the solution u(f.t) and also the inf-sup infimizer at the {i n sample points; and the 
(lngh-dimensional) supremizing space is taken to be X. Stability is ensured and rigorous (sharp) error 
bounds are obtained though technically the bounds are only asymptotic due to the approximation of 
the inf-sup paiameter; and, despite the presence of X, the on-line complexity remains independent of 
the dimension of X as in Section 3.3, we exploit affine parameter dependence and linear superposition 
to precompute the necessary inversions. In a second suite of much simpler and more general approaches 
(see [17] for one example in the symmetric case), we exploit minimum-residual or Petrov- Galerkin pro- 
jections with infimizer -supremizer enriched, but still very low-dimensional, infimizing and supremizing 
spaces. Plausible but not yet completely rigorous arguments, and empirical evidence, suggest that 
stability is ensured and rigorous asymptotic (and sharp) error bounds are obtained. 

In [17] we focus entirely on Method I a posteriori error estimator procedures; but Method II tech- 
niques axe also appiopriate. In particular, Method II approaches do not require accurate estimation of 
the inf-sup paiametei, we thus need be concerned only with stability in designing our reduced— basis 
spaces. 


5.3.2 Parabolic Partial Differential Equations 

The next extension consideied is the treatment of parabolic partial differential equations of the form 
v \ /*) ~ a (' u i v \ typical examples are time-dependent problems such as unsteady heat conduction 
— the “heat” or “diffusion” equation. The essential new ingredient is the presence of the time variable 
t . 

The reduced-basis approximation and error estimator procedures are similar to those for noncompli- 
ant nonsym metric problems, except that we now include the time variable as an additional parameter. 
Thus, as in certain other time-domain model-order-reduction methods [4, 25], the basis functions are 
snapshots of the solution at selected time instants; however, in our case, we construct an ensemble of 
such series corresponding to different points in the non-time parameter domain T>. For rapid conver- 
gence of the output approximation, the solutions to an adjoint problem — which evolves backward in 
time — must also be included in the reduced-basis space. 

For the temporal discretization method, many possible choices are available. The most appropriate 
method — although not the only choice — is the discontinuous Galerkin method [13]. The varia- 
tional origin of the discontinuous Galerkin approach leads naturally to rigorous output bounds for 
Method I a posteriori error estimators; the Method II approach is also directly applicable. Under our 
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affine assumption, off- -line /on-line decompositions can be readily crafted; the complexity of the on-line 
stage (calculation of the output predictor and associated bound gap) is, as before, independent of the 
dimension of X, 

5.3.3 Locally Non- Affine Parameter Dependence 

An important restriction of our methods is the assumption of affine parameter dependence. Although 
many property, boundary condition, load, and even geometry variations can indeed be expressed in the 
required form (2) for reasonably small Q, there are many problems — for example, general boundary 
shape variations — which do not admit such a representation. One simple approach to the treatment 
of this more difficult class of non- affine problems is ( i ) in the off -line stage, store the ( n = u(/x n ), and 
(u) in the on-line stage, directly evaluate the reduced— basis stiffness matrix as /x). Unfortu- 

nately, the operation count (respectively, storage) for the on-line stage will now scale as Q(N 2 dim(X)) 
(respectively, 0{N dim(X)), where dim(X) is the dimension of the truth (very fine) finite element 
approximation space: the resulting method may no longer be competitive with advanced iterative tech- 
niques; and, in any event, “real-time” response may be compromised. 

We prefer an approach which is slightly less general but potentially much more efficient. In partic- 
ular, we note that in many cases — for example, boundary geometry modification — the non- affine 
parametric dependence can be restricted to a small subdomain of ft, ft//. We can then express our 
bilinear form a as an affine/ non -affine sum, 

a(w,v;fi) = a I (w,v\(i) + a n (u>, v\ /x). (59) 

Here a/, defined over ft/ — the majority of the domain — is affinely dependent on /x; and a//, defined 
over ft// — a small portion of the domain — is not affinely dependent on /x. It immediately follows that 
the reduced -basis stiffness matrix can be expressed as the sum of two stiffness matrices corresponding 
to contributions from a/ and a u respectively; that the stiffness matrix associated with a/ admits the 
usual on- line/ off-line decomposition described in Section 3.3; and that the stiffness matrix associated 
with an requires storage (and inner product evaluation) only of (C* restricted to ft//). The 

non-affine contribution to the on-line computational complexity thus scales only as 0(N 2 dim(A"|n /J )), 
where dim(X|n i; ) refers (in practice) to the number of finite-element nodes located within ft// — often 
extremely small. We thus recover a method that is (almost) independent of dim(X), though clearly the 
on-line code will be more complicated than in the purely affine case. 

In the above we focus on approximation. As regards a posteriori error estimation, the non-affine 
dependence of a (even locally) precludes the precomputation and linear superposition strategy required 
by Method I (unless domain decomposition concepts are exploited [12]); however, Method II directly 
extends to the locally non-affine case. 
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